Hartree— Fock Approximation for Inverse Many— Body Problems 
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A new method is presented to reconstruct the potential of a quantum mechanical many-body 
system from observational data, combining a nonparametric Bayesian approach with a Hartree- 
Fock approximation. A priori information is implemented as a stochastic process, defined on the 
space of potentials. The method is computationally feasible and provides a general framework to 
treat inverse problems for quantum mechanical many-body systems. 
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The reconstruction of inter-particle forces from obser- 
vational data is of key importance for any application 
of quantum mechanics to real world systems. Such in- 
verse problems have been studied intensively in inverse 
scattering theory and in inverse spectral theory for one- 
body systems in one and, later, in three dimensions jl],^| . 
In this Paper we now outline a method, designed to deal 
with inverse problems for many-body systems. 

Being the mathematical counterpart of induction prob- 
lems in philosophy, inverse problems appear quite natu- 
rally in science when justification of a physical law has to 
be based on a finite number of observations. Such prob- 
lems are notoriously ill-posed in the sense of Hadamard 
In that case it is well known that additional a pri- 
ori information is required to obtain a stable and unique 
solution. Referring to a Bayesian framework j^-RL we 
implement a priori information in form of stochastic pro- 
cesses over potentials ||] • While a standard procedure is 
to fit parameterized potentials to the data, we will es- 
pecially be interested in less restrictive, nonparametric 
approaches. As, up to now, calculating an exact solution 
of inverse many-body problems is not feasible, we will 
treat the problem in an 'inverse Hartree-Fock approxi- 
mation' (IHFA). 

A main advantage of Bayesian methods is their flexi- 
bility. They can easily be adapted to different learning 
situations and have therefore been applied to a variety 
of empirical learning problems, including classification, 
regression, density estimation ||[l0|], and, recently, to 
quantum statistics |8). In particular, within a Bayesian 
approach it is straightforward to deal with measurements 
of arbitrary quantum mechanical observables, to include 
classical noise processes, and to implement a priori in- 
formation explicitly in terms of the potential. Compu- 
tationally, on the other hand, working with stochastic 
processes, or discretized versions thereof, is much more 
demanding than, for example, fitting parameters. This 
holds especially for applications to quantum mechanics 
where one can not take full advantage of the convenient 
analytical features of Gaussian processes. Due to increas- 
ing computational resources, however, the corresponding 
learning algorithms become now numerically feasible. 

We will consider many-fermion systems with Hamilto- 



nians, H = T + V, consisting of a one-body part T, e.g., 
T = — (1/2to)A (with Laplacian A, mass m,fi = 1), and 
a two-body potential V. To write such Hamiltonians in 
second quantization, we introduce creation and annihila- 
tion operators a) a , a a corresponding to a complete single 
particle basis \ip a >, i.e., at|0> = | tp a > and a a |0> 
= 0. As we are dealing with fermions, we have to re- 
quire the usual anticommutation relations a Q a 7 + ata a 
= < <p a | tp 1 >, a^a^ + a* at =0, a Q a 7 + a 7 a Q =0, and we 
can write 



H 



a/3 
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with T a _p = < (p a \T\ipf3 > and antisymmetrized matrix 
elements Vap-ys — <(p a ( Pp\V\(p'y(ps>- We will consider 
two-body potentials which are local and depend only on 
the distance between the particles x = \x\ — ara|, i.e., 
V XlX2X ' lX ' 2 = v(\xi - x 2 \)(S(x 1 - x[)S(x 2 - x' 2 ) - 8{xx - 
x' 2 )S(x2 — x[)). For unknown function v(x), our aim will 
be to reconstruct this function from observational data. 

To obtain information about the potential, the system 
has to be prepared in a state depending on v. Such a 
state can be a stationary statistical state, e.g. a canonical 
ensemble, or a time-dependent state evolving according 
to the Hamiltonian of the system. In the following we 
will study many-body systems being prepared in their 
ground state. The (normalized) TV-particle ground state 
wave function ipo depends on v and is antisymmetrized 
for fermions. As observational data we choose n simul- 
taneous measurements of the coordinates of the N parti- 
cles, the corresponding observable being the coordinate 
operator x. The ith measurement results hereby in a vec- 
tor Xi , consisting of N components Xij , each representing 
a single particle coordinate. Introducing the Slater de- 
terminant | Xi > = | aft i, ■ • ■ , Xi tll >, made of orthonormal 
single particle orbitals | Xjj >, the probability density of 
measuring the coordinate vector x\ given v is, according 
to the axioms of quantum mechanics, 



p(Xi\x,v) =<1po\Xi><Xi\lp >= \ip (x l 



■ ,Xi,N)\ 2 , 



(2) 



which, when regarded as function of v, for fixed Xi, is also 
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called the likelihood of v. In contrast to an ideal measure- 
ment of a classical system, the state of a quantum system 
is typically changed by the measurement process. In par- 
ticular, its state is projected in the space of eigenfunc- 
tions of the measured observable with eigenvalue equal 
to the measurement result. Hence, if we want to deal 
with independent, identically distributed data, the sys- 
tem must be prepared in the same state before each mea- 
surement. Under such conditions the total likelihood fac- 
torizes 



p(xi, ■ ■ ■ ,x n \x,v) = Y\_p(xi\x,v). 



(3) 



In maximum likelihood approximation, a potential is 
reconstructed by selecting a space of parameterized po- 
tentials v(x, £) and maximizing the total likelihood (|^) 
with respect to the parameters £, i.e., v*(x) = v(x,£*) 
with £* = argmaXf p(x±, ■ ■ ■ , x n \x, «(£))• This, however, 
only yields a unique solution if the parameterized space 
is small enough. Otherwise, additional constraints have 
to be included to determine a potential uniquely. In a 
Bayesian framework those constraints are provided by 
additional a priori information and implemented by se- 
lecting a prior density p(v), interpreted as probability 
density before having received data. The posterior den- 
sity p(v\D), being the probability density of v after hav- 
ing received data D, is then obtained according to Bayes' 
rule, 



p(v\D) = 



P{ v )Y\iP{xi\x,v) 
p(xi, - ■ ■ ,x n \x) 



(4) 



Finally, the predictive density is obtained from the pos- 
terior density as posterior expectation of the likelihood 
p(x\x 1 D) = Jdv p(x\x,v) p(v\D). Within a nonparamet- 
ric approach, where function values v(x), and not pa- 
rameters, represent the fundamental variables, p{v) and 
p(v\D) are a stochastic processes and the integration over 
v is a functional integration over functions v(x). Such in- 
tegrations may be approximated by Monte-Carlo meth- 
ods, or, as we will do in the following, be treated in max- 
imum a posteriori approximation, a variant of the saddle 
point method. In that case one assumes that the main 
contribution to the integral comes from the potential with 
maximal posterior, i.e., p(x\x 7 D) w p(x\x,v*) where v* 
= &Tgmax v p(v\D). So we are left with maximizing the 
posterior (||), which we will do by setting the functional 
derivative of the posterior with respect to v(x) to 
zero (for x > 0). Hence, introducing the notation 8 V < X \ 
= 8 / 8v{x), we have to solve the stationarity equation 



= 8 v(x) \np(v) + ^2 S v(x) lnp(xi\x, v) 



(5) 



but much more general, we can also consider mixtures of 
Gaussian processes [flT|], for which p(v) = ^2f.p(k)p(v\k) 
with Gaussian components 



p(v\k) 



dct • 



2tt 



-±<v-v k \K k \v-v k > 



(6) 



positive (semi-) definite covariance K^ -1 , and regression 
function ffc(x), playing the role of a reference potential. 
A typical choice for the inverse covariance is the negative 
Laplacian multiplied with a 'regularization parameter' A, 
i.e., K/j = — AA, which favors smooth potentials. Higher 
order differential operators may be included in Kfc, as it 
is often done in regression problems to get differentiable 
regression functions 12 Also useful are integral op- 



erators, for example, to enforce approximate periodicity 
of v ||. The functional derivative of a Gaussian mixture 
prior with respect to v is easily found as 



S v (x) lnp(v) = -K (v - v ), 



(7) 



where K = J2k p(M v ) K fe and v o = K 1 J2k Pi k \ v ) K fe«fe 
with p(k\v) — p(v\k)p(k) / p(v) . 

To calculate the functional derivative of the likeli- 
hood @, 5 v ( x )p(Xi\x,v) =< 5 v ( x )4> \xi >< Xi\ip > + 
<ipo\xi><Xi\5 v ( x )%l)Q>, we need 8 v ( x )ipo- It is straight- 
forward to show, by taking the functional derivative of 
HtjjQ — EqiPq, that, for nondegenerate ground state, a 
complete basis of eigcnstates ip~ f with energies E^, and 
requiring <ipo\8 v ( x )ij)o> — to fix norm and phase, 



\S v ( x )ipo>= ^2 



1 
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EQ — Er 



■\ip 1 ><ip 1 \6 v{x) H\<ipo>. (8) 



Furthermore, from 8 V / X \v(\xi — a^l) = 8(x — \x% — a^l) 
directly follows 

$v(x)H = - ^a Xl (a xl _ x a xx - x + a Xl+x a Xl+x )a Xl . (9) 



Typically, a direct solution of the many-body equation 
d|) is not feasible. To get a solvable problem we treat 
the many-body system in Hartree-Fock approximation 
@-@ (For non-hermitian H see |17| ), Thus, as first 
step of an IHFA, we approximate the many-body Hamil- 
tonian H by a one-body Hamiltonian h defined self- 
consistently by 



A' 



T xx > + x(fk I V I x'cpk > , 



(10) 



fe=i 



ifk being the TV-lowest orthonormalized eigenstates of h, 
i.e., 



The most commonly used prior processes are Gaus- 
sian. Being technically only slightly more complicated 



ei < £fc < £at- 



(11) 
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The corresponding normalized TV-particle Hartree-Fock 
ground state <&o = \<pi, • • • } <Pn > is the Slater determi- 
nant made of the TV-lowest single particle orbitals ipk- 
The Hartree-Fock likelihood, replacing (En, becomes, 



,(r+l). 



pn F (xi\x,v) = <$ I Xi ><Xi | $ >- 



(12) 



To find the functional derivative 5 v r x \pHF(xi \x, v) — 

<6 v ( x <)$o\Xi><Xi\$o> + <$0 \xi><Xi |^( x )*0>, 

we define the overlap matrix Bi with matrix elements 
Bki.i = <Pk(zi,l) = < x u I ^fe > and expand <Xj | $ > = 



det Bi = Mu-,%Bu\ii m terms of its cofactors Mu,i — 
(B~ 1 )jfc det B t . Applying the product rule yields 



iV 



$v(x)< x i \$Q>=) j Mkl-i S v(x ) Pk(Xi,l). 



(13) 



Again, we proceed by taking the functional derivative 
of Eq. ( |TT| ) and obtain after standard manipulations (for 
nondegenerate Ck and < ipk\S v (x) l Pk > — 0), 

\S v (x)Vk >= ^2 — >< l Pi\ S v(x)h\(pk > , (14) 



and, following from Eq. ( |l0|) 

N 



5 v (x)hx>x" = ( <x' (fj 1 5 V ( X )V | x" <pj > 

+ <x'8 v ( x )<pj \ V\x" ipj> 

+ <x'cpj | V \ x" 6 v ( x )(fj > y (15) 

Finally, the functional derivative of the orbitals is ob- 
tained by inserting Eq. ( |l5j ) into Eq. ( O ) 



<^(x)¥>fc(z') = ~I ' X! ( <¥W I ^ I <Pk<Pj > 



N 



+ <<Pi8 v ( x )<Pj | V | yjfe^j > 
+ <<W \ V\<pkS v ( x )<pj> 



(16) 



which can quite effectively be solved by iteration, start- 
ing for example with initial guess S v Mifj(x') = 0. Be- 
cause 5 v ( x )<pk(x') depends on two coordinates x and x', 
Eq. (|l6|), being the central equation of the IHFA, has 
the dimension of a two-body equation for the lowest 
TV orbitals. Introducing, analogously to 73j, the matrix 
Ai(x) with Aki-i(x) = <x ii i\6 v ( x )(p k > = Sv(_x)¥k(xi,l) 
(for x > 0, 1 < k, I < TV, 1 < i < n), it is straightforward 
to check that the functional derivative of a likelihood 
term is given by 

S <x) ln PliF (xi\x,v) = 2Re [Tr^A^x))] . (17) 

The stationarity equation (^|) can now be solved by iter- 
ation, for example, 



r)[v - 



,0). 



K 



5 x lnp B _ F (x l \x 1 v ir) )]. 



(18) 



choosing a positive step width rj and starting from an 
initial guess . 

In conclusion, reconstructing a potential from data by 
IHFA is based on the definition of a prior process for v 
and requires the iterative solution of 

1. the stationarity equation for the potential (|5|), need- 
ing as input for each iteration step (18) 

2. the functional derivatives of the likelihoods ( [HI) , 
obtained by solving the (two-body-like) equation ( |16| ) 
for given 

3. single particle orbitals, defined in ( pr|) as solutions 
of the direct (one-body) Hartree-Fock equation (|TT|). 

We tested the numerical feasibility of the IHFA for a 
Hamiltonian H = -5^ A + V x + V (with m = 1(T 3 ), 
including a local one-body potential with Vi(z,z') = 
S(z — z')az 2 (and a = 10 -5 ) to break translational sym- 
metry. We want to reconstruct the unknown local two- 
body potential V from empirical data. To be able to 
sample data from the 'true' many-body likelihood (|^) 
and to check the quality of the IHFA for a given 'true' 
potential Vtruo, we have to solve the corresponding many- 
body problem. Therefore, we have chosen a system of 
two particles on a one-dimensional grid (with 21 points) 
for which the true ground state can be calculated by 
diagonalizing H numerically. We want to stress that 
application of the IHFA to systems with TV > 2 par- 
ticles is straightforward and only requires to solve Eq. 
( p6| ) for TV instead for two orbitals. Hence, selecting 
a 'true' local two-body potential Vt ruo with Vt TU c{x) = 
b/{l + e- 2 ^ x - L l 2 ^ L ) (and b = 100, 7 = 10, L = 21), we 
were able to sample 100 data points from the correspond- 
ing 'true' probability density (0). The true likelihood as 
function of inter-particle distances x and the empirical 
density of distances p em p(x) = (1/n) Yh=i S x-\xi,j-x it2 \ 
obtained from the training data are shown in Fig. |l|. 

We used a nonparametric approach for v, combined 
with a Gaussian smoothness prior with inverse covariance 
K =A(I-A)/2 (identity I, A= lO" 3 / 2 ), and a reference 
potential vq(x) of the form of ftrue, but with 7=1 (so 
it becomes nearly linear in the interval [1,20]). Further- 
more, we have set all potentials to zero at the origin and 
constant beyond the right boundary. The reconstructed 
potential wihfa has then been obtained by iterating ac- 
cording to Eq. ( |l8| ) and solving Eqs. (^lj) and ( |l6| ) within 
each iteration step. The resulting IHFA likelihood pihfa 
= p(x\x, ^ihfa) indeed fits well the true likelihood ptruo 
= p(x\x,Vi; Iue ) (see Fig. |l|). In particular, pihfa is over 
the whole range an improvement of the reference likeli- 
hood po = p(x\x,Vq). The situation is more complex for 
potentials (see Fig. |^). On the basis of 100 data points, 
the true potential is only well approximated at medium 
inter-particle distances. For large and small distances, 
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on the other hand, the IHFA solution is still dominated 
by the reference potential of the prior process. This ef- 
fect is a consequence of the lack of empirical data in those 
regions (see Fig. [j]): The probability to find particles at 
large distances is small, because the true potential has 
its maximum at large distances. Also, measuring small 
distances is unlikely, because antisymmetry forbids two 
fermions to be at the same place. In such low data re- 
gions one must therefore rely on a priori information. 

We are grateful to A. Weiguny for stimulating discus- 
sions. 
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FIG. 1. Empirical density p e m P for n = 100, true likelihood 
Ptruo, reference likelihood po, reconstructed likelihood pihfa 
as function of inter-particle distance. 
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Inter-particle distance 

FIG. 2. True potential utrus, reference potential vo, and re- 
constructed IHFA potential uihfa as function of inter-particle 
distance. 



4 



